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Land use change in the form of urbanization is a direct driver affecting 
the provision of ecosystem services from forests. To better understand 
this driver, we modeled the effects of urbanization on three regulating 
and provisioning ecosystem services in two disparate watersheds in 
Florida, USA. The study integrated available geospatial and plot-level 
forest inventory data to assess future changes in carbon storage, timber 
volume and water yield during a period of 57 years. A 2003-2060 
urbanization and land use change scenario was developed using land cover 
data and a population distribution model. The Integrated Valuation and 
Ecosystem Services Tradeoffs model was then used to quantify changes in 
ecosystem services. Carbon storage was reduced by 16% and 26% in the 
urbanized 2060 scenario in both the rural Lower Suwannee and urban 
Pensacola Bay watersheds, respectively. Timber volume was reduced by 11% 
in the Lower Suwannee and 21% in the Pensacola Bay watershed. Water 
yield, however, increased in both watersheds by 4%. Specific sub¬ 
watersheds that were most susceptible to urbanization were identified and 
mapped and ecosystem service interactions, or trade-offs and synergies, 
are discussed. Findings reveal how urbanization drives the spatio- 
temporal dynamics of ecosystem services and their trade-offs. This study 
provides policy makers and planners an approach to better develop 
integrated modeling scenarios as well as designing mapping and monitoring 
protocols for land use change and ecosystem service assessments. © 2016 
Elsevier Ltd. All rights reserved. 1. Introduction The world's population 
is rapidly increasing (US Census Bureau, 2013) and the world's population 
is projected to be 9 billion by 2050. This population increase will lead 
to land use and cover (LULC) changes and alter the provision of ecosystem 
services such as food, timber, and clean water (European Commission, 

2009; Hoyer and Chang, 2014; Millennium Ecosystem Assessment, 2003; 
Vitousek et al., 1997). Urbanization, the development of new urban areas 
from non-urban lands, is a key anthropogenic driver affecting ecosystems 
(Bengston et al., 2005; He et al., 2016; Yue et al., 2013; Zhang et al., 
2012). The use of the driver concept in influenc- * Corresponding author. 
Current address: Functional and Ecosystem Ecology UnitBiology Program, 
Faculty of Natural Sciences and Mathematics, Universidad del Rosario, Kr 
26 No 63B-48, Bogota D.C., Colombia. E-mail addresses: 
soniadelphin@gmail.com (S. Delphin), fjescobe@gmail.com, 
franciscoj.escobedo@urosario.edu.co (F.J. Escobedo), aamr@ufl.edu (A. 
Abd-Elrahman), wcropper@ufl.edu (W.P. Cropper). 1 Current address: World 
Wildlife Fund, Edificio Opa Rudy-4to piso, Avda. Argana~ casi Avda. 

Peron, Asuncion, Paraguay, ing land use change and ecosystem services, 
has become a common approach for studies analyzing policies andmanagement 
objectives and associated effects on human well-being (Hoyer and Chang, 
2014;Yue et al., 2013) .By analyzing ecosystemservicedriversusing 



integrated modeling scenarios, as well as measured and available data, we 
can better understand ecosystem function and changes in their services. 

We can also provide information for planning and policy formulation in 
order to minimize potential negative impacts to human well-being 
(Millennium Ecosystem Assessment, 2003) . Forests are important for carbon 
storage and this function is referred to as a regulating ecosystem 
service as it regulates global climate (Millennium Ecosystem Assessment, 
2003). Forest carbon storage is commonly measured according to four 
pools: aboveground biomass, belowground biomass, litter or dead biomass, 
and soil (Ashton et al., 2012) . Indeed, forests are critical to the 
global carbon cycle as they comprise a large stock of carbon relative to 
other ecosystems; thus conserving peri-urban forested areas from 
urbanization is a priority (Gorte, 2009; Harmon, 2001; Jandl et al., 

2007; Yonavjak et al., 2011). Urbanization of forests in the southern 
U.S. has for example resulted in about 0.21 Petagrams (Pg) of carbon 
emissions from 1945 to 2007 (Zhang et al., 2012). Zhang and 
http://dx.doi.Org/10.1016/j.landusepol.2016.02.006 0264-8377/© 2016 
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the principal drivers in the decline of natural forest timber areas. 

Urban areas can store carbon in the form of planted and remnant trees; 
but, 70-100 years are needed before pre-urbanization carbon stocks are 
reached (Zhang et al., 2012). Furthermore, the recovery of carbon stocks 
to pre-urbanization levels is uncertain and depends on the specific 
ecosystem, its structure and management of the existing urban forest. As 
such, timber supplies are significantly affected by urbanization. Timber 
products are considered a provisioning ecosystem service that provides 
direct benefits to people (Boyd and Banzhaf, 2007; Millennium Ecosystem 
Assessment, 2003) . The benefits of these provisioning ecosystem services 
(e.g., timber, food, nontimber forest products, water supply) are widely 
recognized and are more easily valued than most other ecosystem services 
because of the existence of markets (De Groot et al., 2010) . For example, 
forests in the Southern U.S. produce more than 60% of the nation's timber 
products and substantially contribute to the region's economy (Prestemon 
and Abt, 2002a). Bystriakova et al. (2005) report that timber harvesting 

at the global scale has increased by 60% in the last four decades and the 
increasing trend is expected to continue, but at a slower rate. Further, 
Cubbage et al. (1995) found that urbanization reduced timber supply, 

especially of softwoods and similarly Hodges et al. (1998) has documented 

how forest-timber areas in Louisiana, U.S. have been reduced by urban 
development. Water is another valuable forest regulating and provisioning 
ecosystem service (Kreye et al., 2014). Land use changes can however 
alter the quantity and quality of water (Foley et al., 2005; Sun et al., 
2005); therefore, determining the effect of forest LULC changes on water 
yield and quality is an important consideration. Water yield is the sum 
of annual precipitation that does not evaporate from soil and water or 
transpire from vegetation (Mendoza et al., 2011). Since urban areas 
increase runoff, as a result of decreased interception, 

evapotranspiration and infiltration, they can decrease water quality and 
increase overall water yield (Arnold et al., 1987; Hanson et al., 2010). 
Indeed several studies have determined the relationship between the 
hydrological cycle, forests (Bosch and Hewlett, 1982; Sun et al., 2005), 
and urbanization. Hollis (1977) reports that median flow increased by 
approximately 40% post- urbanization in the Canon's Brook catchment in 
England. Arnold et al. (1987) also modeled the effect of urbanization on 

water yield in a rural Texas U.S. watershed that urbanized over 70 years 
to about 77% urban land cover and found that annual surface runoff would 
increase by about 10%. Other studies in Florida, U.S. have analyzed the 
effect of urbanization on water quality (Cooley and Martin, 1979; Frick, 
1998) and water supply (Boggess, 1968). However, we found no other 



relevant studies on the effects of urbanization on forest water yield or 
timber production in the State of Florida. Geographic Information Systems 
(GIS) models have been used to assess the potential effects of different 
drivers on the provision of ecosystem services in subtropical forests 
(Cademus et al., 2014; Metzger et al., 2006). Timilsina et al. (2011), 

for example, used geospatial models to map forest carbon storage hotspots 
in Florida U.S. and identified plot-level ecological and anthropogenic 
drivers ofthis regulating ecosystem service, based on the Millennium 
Ecosystem Assessment (2003) classification. Cademus et al. (2014) also 

used geospatial and forest inventory data to assess interaction or the 
trade-offs (win-lose situations), synergies (win-win situations) and 
drivers in pine dominated forests in Florida. Models such as the 
Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) 
developed by the Natural Capital Project (http://www. 

naturalcapitalproject.org/), have also used available geospatial data to 
estimate and map the provision of ecosystem services (Tallis et al., 

2011). The model estimates the provision of ecosystem services in a 
current land use/cover and in different scenarios that incorporate 
proposed policies, management goals and disturbance regimes. Because of 
this, the InVEST model has been recently used to quantify and value 
ecosystem services and thus facilitate their integration into policy 
making (Daily et al., 2009; Hoyer and Chang, 2014; Tallis et al., 2011). 
Daily et al. (2009) and Nelson et al. (2009) have applied InVEST in 
different areas of the US, and He et al. (2016) in China, to assess the 

effects of alternative land use and management decisions on ecosystem 
services. Hoyer and Chang (2014) estimated the provision of freshwater 
ecosystem services for two basins in Oregon using several urbanization 
and climate change scenarios. Delphin et al. (2013) and Bottalico et al. 

(2016) also analyzed the effect of a driver (i.e., hurricanes) and forest 
management scenarios, respectively, on two ecosystem services,timber and 
carbon stocks, using the InVEST and other geospatial models. Approaches 
such as those of Delphin et al. (2013), Bottalico et al. (2016), He et 

al. (2016); and Hoyer and Chang's (2014) could therefore be adapted to 

assess both the spatial and temporal effects of urbanization on multiple 
ecosystem services and their interactions as well as to better formulate 
land use policies. Given the lack of studies on the effects of ecosystem 
service drivers on subtropical forest ecosystem service provision, the 
aim of our study was to assess the effects of urbanization on a bundle of 
forest ecosystem services and their trade-offs: carbon storage as an 
indicator of climate regulation, water yield as an indicator of flooding 
regulation, and timber volume, as an indicator of a provisioning service. 
Our specific objectives were to quantify both the temporal and spatial 
effects of urbanization on these key ecosystem services from 2003 to 2060 
using: (1) the InVEST model, (2) plot-level United States Department of 
Agriculture (USDA) Forest Inventory and Analysis (FIA) data, and (3) a 
regional forecasted urbanization land cover change scenario for 2003-2060 
based on a population distribution spatial model. We used the rural Lower 
Suwannee River and urbanized Pensacola Bay watersheds in northwestern 
Florida U.S. as our two study areas. These two representative, yet 
disparate in terms of rural versus urban land uses, study watersheds were 
selected based on: the proportion of forest and urban land covers, FIA 
plot density, and the presence of impaired water bodies (i.e., do not 
meet water quality standards) according to Florida's Statewide 
Comprehensive Verified List of Impaired waters 

(http://www.dep.state.fi.us/ water/watersheds/assessment/a-lists.htm). 
This integrated modeling approach using two subtropical watersheds should 
contribute to policy-relevant land management and planning assessments. 

It can also be used to better communicate the temporal and spatial 
effects and trade-offs of urban development and other land use changes on 
ecosystem services in other types of forested ecosystems in the world. 2. 



Material and methods 2.1. Study areas The rural Lower Suwannee River 
watershed in western Florida has 45% forest, 6% urban, and 25% 
agricultural land covers while the urbanized Pensacola Bay watershed in 
northwestern Florida has approximately 20% urban, 39% forest, and 0.1% 
agriculture land covers (Fig. 1) . The Lower Suwannee watershed consists 
of 63 subwatersheds in an area of 408,828 ha. and is one of five 
watersheds that comprise the Suwannee River watershed (Katz et al., 

1997). The Suwannee river is the second largest river in Florida in terms 
of average discharge (Light et al., 2002). The Pensacola Bay watershed 
comprises 15 sub-watersheds in an area of 140,825 ha. and is part of the 
Pensacola Bay system that is composed by five watersheds (Schwenning et 
al., 2007) . As previously explained, the two 190 S. Delphin et al. / Land 
Use Policy 54 (2016) 188-199 Fig. 1. The Lower Suwannee River and 
Pensacola Bay watershed study areas in Florida U.S. watersheds were 
selected as they have been prioritized by the State of Florida for water 
quality protection. 2.2. Geospatial and field data 2.2.1. Florida 
vegetation and land cover 2003 We used the 43 vegetation and land cover 
classes from the 2003 Florida Vegetation and Land cover (FVLC) map 
developed by the Florida Fish and Wildlife Conservation Commission (FWC) 
based on Landsat Enhanced Thematic Mapper satellite imagery with a 30 
meter spatial resolution (Stys et al., 2004) . Carbon stores and timber 
volume were estimated for the 11 forest cover classes identified by FVLC 
data. Given that the scope of this analysis is forest land cover changes, 
forests or trees and other associated vegetation within already urbanized 
land covers, were not spatially or statistically analyzed. However, 
existing studies were used to explore carbon and water yield dynamics 
within urban forests and land covers. At the time of this study the 2003 
FVLC map was the most detailed classification for the state of Florida. 
2.2.2. 2060-Land use and land cover map The forecasted urbanization 2060- 
LULC map was developed from the 2003-Florida Vegetation and Land Use map 
and a 2060- Florida Population Distribution model (Zwick and Carr, 2006). 
This model is based on population projections from the Bureau of Economic 
and Business Research (BEBR) and a GIS land use suitability analysis 
(Zwick and Carr, 2006). The data and models account for periods of 
increased economic activity and rapid urbanization (i.e., 2003-2008), 
economic downturn and decline in urbanization rate (i.e., 2008-2013), and 
planning projections beyond 2040 that are regularly used for urban and 
forest policy formulation in Florida (Wear and Greis, 2002) . Three main 
criteria considered in the land use model were: (1) suitability to 
support a specific land use category, (2) preference to reflect community 
values, and (3) potential conflict with conservation (i.e., forests), 
urban and agricultural lands. Based on BEBR population projections, the 
model calculates the human population and land area that need to be 
allocated; it then identifies and maps suitable land use polygons for 
these new urban areas in Florida. Overall, polygons that are close to 
existing urban areas, higher road densities, and that lack wetlands and 
conservation lands were the ones most heavily weighted for land use 
change. Specific details, methods and assumptions can be found in Zwick 
and Carr (2006) . The 2060-Population Distribution map contains the 
existing and projected urban areas, this result was then combined with 
the 2003-Florida Vegetation map to forecast and produce the 2060 
urbanization scenario-land cover map. 2.2.3. Forest inventory and 
analysis program data We used available USDA Forest Service Inventory and 
Analysis (FIA) program data to estimate forest carbon storage and timber 
volume as the 2003 current baseline for our urbanization scenario. The 
FIA provides tree and plot-level data for all forest lands in all the U.S 
including forest extent, characteristics, condition, and tenures (Jenkins 
et al., 2003; Woudenberg et al., 2010) . We specif- S. Delphin et al. / 
Land Use Policy 54 (2016) 188-199 191 ically used the FIA data from cycle 
8 (2002-2007) since it accounts for our analysis period. To ensure 



private landowners' privacy and plot integrity, geographic coordinates 
(latitude and longitude) of the plots are not provided because of a 
privacy provision enacted by the U.S, Congress in the Food Security Act 
of 1985. Therefore, plots are located within a buffer of approximately 
0.8 to 1.6 km from their true centroid. Alternative approaches based on 
other FIA attributes, can however be applied to better approximate true 
plot centroids. For example, only the 1-10% of plots are misclassified 
when analyzing for ecological subsection data (Sabor et al., 2007). Thus, 
accounting for the ecological subsection strata, FIAsurvey unit and 
County attributes can better help specify true plot centroids. 
Accordingly, we used an approach described in Delphin et al. (2013) that 

uses specific FIA database attributes: FIA Survey Unit, ecological 
subsection, and county, to better approximate the true FIA plot centroid, 
according to unique combination of classes for each forest ecosystem, 
ecological subsection and county. Once forest ecosystems were identified 
for each watershed, the main forest tree species per watershed were 
identified based on vegetation descriptions prepared by Gilbert and Stys 
(2004) that contain the list of forest tree species observed in each 
forest ecosystem. Finally, the forest tree species found in the FIA 
database were considered and the species codes were identified. 2.2.4. 
Forest carbon storage estimates The InVESTCarbonStorage 
andSequestrationmodel wasusedto analyze carbon storage in the 2003- 
baseline land covers. These carbon storage estimates were then projected 
to the 2060-urbanized land cover. The FIA tree and plot-level data 
described in Section 2.2.3 wereusedas inputfor themodel, whichestimates 
carbonstorage in Megagrams (Mg) per hectare (ha)for each land cover. 
Carbon storage analyses were estimated for the differentforesttypes in 
the 2003 Florida Vegetation and Land cover map. The four carbon pools (Mg 
per ha) were aboveground biomass, belowground biomass, soil organic 
matter and dead biomass (Appendix A; Timilsina et al., 2011). The FIA 
provided above and belowground carbon stores at the individual tree- 
level; therefore to obtain plot-level data, per-tree values were 
multiplied by an FIA-defined expansion factor to scale trees on a plot to 
their peracre equivalent (Woudenberg et al., 2010). Accordingly, carbon 
values per tree were summed by plot using species code (SPCD), FIA survey 
unit (UNITCD), ecological subsection (ECOSUBCD), county (COUNTYCD), 
condition status (COND STATUS CD), and stand age (STDAGE). Only plots 
with the condition status defined as forest land and with a STDAGE 
ranging from 1 to 72 years were considered in the analysis. We note that 
STDAGE is not the overall age of a stand but is an FIA indicator of the 
average live age of trees in a plot that are not overtopped in the 
predominant size-class for a given condition. Please see Appendix A for a 
description of carbon pools and Woudenberg et al. (2010) for a definition 

of FIA codes and variables. Timber harvest will reduce the amount of 
carbon stored in forests, however, the specific impact on landscape level 
carbon storage will depend on the type of silvicultural approach used 
(e.g., less frequent harvesting) (Schwenk et al., 2012) . Because we were 
not able to account for this, no specific timber harvest activity effects 
were considered in the analysis. The FIA provides plot-level dead biomass 
and soil carbon estimates in tons per acre, which were then converted to 
megagrams per hectare. The carbon pools for the urbanized 2060 scenario 
were also estimated using an assumption that carbon pools in 2060 will be 
5 percent less than 2003 pools (Huggett et al., 2013). There is 
significant uncertainty in future forest carbon projections, and so based 
on USDA Forest Service forecasts, this 5% reduction in annual forest 
growth for the southeastern region of the US was incorporated into our 
estimates (Heath and Smith, 2000; Huggett et al., 2013; Moore et al., 

2002) . This assumes that carbon stock estimates are a simple linear 
factor and thus they can easily be adjusted. The InVEST carbon storage 
model's intermediate output was the sum of carbon stocks in the four 



pools: carbon aboveground (all living plant material above soil), 
belowground (living root systems below soil), dead (litter and standing 
dead biomass), and soil. Final model output was expressed in total Mg of 
carbon per pixel that were then converted to Mg per hectare. Data was 
classified using the ArcGIS Natural Breaks method (ESRI, 2011) to display 
high, moderate and medium range of change (2003-2060) for each sub¬ 
watershed. 2.2.5. Timber volume estimates Methods for estimating timber 
production were similar to those for carbon storage. Specifically, timber 
volume and timber in the sawlog portion were calculated as indicators 
ofthe potentialtimber available for harvesting. Plot-level timber volume 
was estimated using FIA survey unit, ecological subsection, and county. 
From the FIA tree-level data, two attributes were analyzed: the net 
volume (volnet) and the net volume in the saw-log portion (volsaw). The 
final results were converted to cubic meters per hectare and volume 
values were assigned to individual pixels in the map. For the 2060 timber 
volume projections, we used the volume values presented in Prestemon and 
Abt (2002b) and the Subregional Timber Supply Model (SRTS) of Abt et al. 
(2000), assuming an annual average increases of 1.7% to determine 2006 
timber volume. Timber data were also classified based on the ArcGIS 
Natural Breaks method (ESRI, 2011)to identify high, moderate and medium 
ranges of change (2003-2060) for each sub-watershed. 2.2.6. Water yield 
estimates Water yield was estimated for a 10 year period using annual 
average precipitation (P) data for 1998-2008, from the PRISM Climate 
Group (PRISM Climate Group and Oregon State University, 2008). For 2060, 

P was estimated using the A2 scenario from Coulson et al.'s (2010) 
projections from the Australian Commonwealth Scientific and Industrial 
Research Organization's Mk3Climate System. The A2 scenario is 
conservative for Florida as it represents lower socioeconomic and urban 
development rates, thus less conversion of forests to urban areas. 

Maximum soil and water table depth required by InVEST were kept constant 
for the analysis period and were obtained from the U.S. General Soil Map 
(State Soil Geographic-STATSG02) developed by the USDA Natural Resources 
Conservation Service (USDA-NRCS) in 2006. Average annual Potential 
Evapotranspiration (PET) for the analysis period corresponded to 1998- 
2008 and was from the U.S. Geological Survey's Florida Integrated Science 
Center and plant available water content (AWC) data for the analysis 
period were from the USDANRCS' General Soil Map-STATSG02. Finally, 
maximum root depth was collected from Canadell et al. (1996) and the 

Zhang Constant (Z), a factor used to characterize the seasonality of 
precipitation in an area was from Tallis et al. (2011) . Using this input 

data, the InVEST Water Purification model estimated the annual average 
water yield for a specific watershed. The model estimated water yield at 
the sub-watershed and watershed level and accounts for surface water and 
not ground water processes (Mendoza et al., 2011) . Thus, there is no 
differentiation between surface, subsurface and base flow water in the 
InVEST model(Tallis et al., 2011) . According to Zektser and Loaiciga 
(1993), groundwater comprises on average 30% of the stream flow in most 
hydrologic systems and in some areas can reach up to 90% (Winter et al., 

1998) . Therefore, in areas where groundwater is influential as in Florida 
US, the InVEST model can produce low water yield values when compared to 
observed stream flow levels (Mendoza et al., 2011). 192 S. Delphin et al. 
/ Land Use Policy 54 (2016) 188-199 The InVEST water yield model is based 
on the Budyko curve and annual average precipitation (Tallis et al., 

2011). The Budyko curve describes the relationship between precipitation 
and actual PET (Mendoza et al., 2011) . This approach assumes that 
evapotranspiration depends on water availability (i.e., precipitation), 
atmospheric demand (i.e., PET; Mendoza et al., 2011; Zhang et al., 2001, 

1999) , and is used to explain the dependence of annual water balance on 
water and energy availabilities (Yang et al., 2011) . The annual water 
yield is therefore estimated based on precipitation and actual 



evapotranspiration. Specifically,themodel analyzes eachpixel considering 
each land cover according to Eq. (1) Yxj = j x 1 - AETxj Pxj x Pxj x 

Axj (1) where Yxj is the annual water yield, AETxj is the annual actual 
evapotranspiration on pixel x with LULCj,and Pxj is the annual 
precipitation in pixel x with LULCj, and Axj is the area of x in LULCj, 
and accordingly AET, P and A were calculated for each LULC type. The 
evapotranspiration partition of the water balance is also calculated. 

This is an approximation of the Budyko curve according to Zhang et al. 
(2001) and is based on PET and plant available water content. The maximum 
difference in AET/P ratio corresponds to forest cover areas where treed 
areas have greater soil water storage (Zhang et al., 2001) and is 
described as: AETxj Pxj = 1 + wxjRxj 1 + wxjRxj + 1/Rxj (2) 

where Rxj is the ratio of potential evapotranspiration to precipitation 
(i.e., Budyko dryness ratio) and wxj is a modified dimensionless ratio of 
plant accessible water storage relative to expected precipitation during 
the year (Mendoza et al., 2011). wxj is given by Eq. (3): wxj = Z AWCx xj 
(3) where AWCx is the volumetric (mm) plant available water content, Z is 
the seasonality factor for rainfall distribution (Mendoza et al., 2011). 
Finally, the ratio of PET to precipitation (Rxj) is calculated by Rxj = 
kj • ETox Pxj = PETxj Pxj (4) where ETox that is the reference 
evapotranspiration from pixel x and kxj is the evapotranspiration 
coefficient associated with LULCj on pixel x.We note that Sanchez-Canales 
et al. (2012) observed that Z was less important than other inputs whereas 
plant evapotranspiration for each LULC was potentially a very sensitive 
parameter in the InVEST model. Water yield (2003-2060) was also 
classified using the ArcGIS Natural Breaks method (ESRI, 2011) to 
identify sub-watersheds with high,moderate andmediumchanges in water 
yield. The water yield model performance was also assessed by comparing 
average annual water yield from the model to average annual stream flow 
data from the Gopher River gage station (site 02323592) located near the 
Lower Suwanee watershed's outlet. Normally, 10 year average values of 
streamflow data shouldbeused for water yield estimates. However, U.S. 
Geological Service (USGS) NationalWater Information System 
(NWIS)(http://waterdata.usgs. gov/nwis/), stream flow data were only 
available for 2000-2008 (Nielsen and Norris, 2007) . 2.3. Statistical 
analysis Statistical differences between the 2003 and the 2060 
urbanization maps at the sub-watershed level were determined using a two 
tailed t-test with two paired samples for means at a 95% confidence 
level. Pearson correlation coefficients were calculated between carbon 
storage, timber volume, water yield, and the amount of forest areas. 
Negative correlation coefficients were also used to determine tradeoffs 
while positive coefficients determined synergies (Cademus et al., 2014; 
Raudsepp-Hearne et al., 2010). All statistical analyses were at the sub¬ 
watershed level using Microsoft Excel. Additionally, three-way trade-offs 
and synergies among carbon storage, timber volume, and water yield 
provision levels were then graphed and assessed using R's (Version 3.0.2; 
http://www.Rproject.org/) scatterplot3d package. 3. Results Findings 
indicate an overall reduction in total carbon storage and timber volume 
and an increase in water yield in both the Pensacola and Lower Suwannee 
river watersheds as a result of urbanization and associated forest cover 
loss. Indeed, these spatiotemporal changes were different in the two 
watersheds given their inherent land use-cover compositions during the 
2003 to 2060 analysis period. As we detail below, ecosystem drivers 
affected all three services, except for timber in the Pensacola 
watershed. Tradeoffinteractions among the three ecosystemservices were 
also evident during the analysis period except for a weak interaction 
between water yield and both carbon storage and timber volume in 2060 in 
the Pensacola Bay watershed. 3.1. Land use and land cover map in 2003 and 
2060 The Pensacola Bay watershed was more urbanized than the Suwannee 
watershedin2003.As expected, PensacolaBay wasmore impacted by 



urbanization in 2060. There was a reduction in forest areas of 21%due 
tourbanization, whereas the Lowe SuwanneeRiver watershed experienced a 
decrease of 12% of forest areas due to the development of new urban 
areas. 3.2. Carbon storage in 2003 and urbanized 2060 The total carbon 
storage from 2003-2060 was reduced by 19% in the Lower Suwannee River 
(Table 1) and this corresponded with a 12% reduction in forest area in 
the urbanized 2060 scenario. There was a reduction in mean carbon stores 
as a result of urbanization across all 63 sub-watersheds. The largest 
difference was observed in the Waccasassa Slough sub-watershed (Id 49, 
Fig. 2), where the mean carbon storage decreased from 183 to 156 Mg/ha 
(15%) as a result of a 14% decrease in forest area and 40% increase in 
urban areas. Fig. 2 shows the range of changes in high, moderate, and low 
change categories. Differences inmeancarbonstorage between all the same 
sub-watersheds in 2003 and 2060 were statistically significant. In the 
Pensacola Bay, the total carbon storage was reduced by 26%during 2003- 
2060 (Table 1) and forest area was reduced by 21% - assuming 5% less 
carbon - in 2060. The mean carbon storage was 148 Mg/ha in 2003 and 137 
Mg/ha in 2060. At the sub-watershed level, the largest difference was 
observed in the Fundy BayouWilliams Creek Frontal sub-watershed (number 
11, Fig. 3), where mean carbon storage dropped from 162 Mg/ha to 144 
Mg/ha; forest area was reduced by 34% and urban area increased by 36%. 
There was a statistically significant reduction of 3-18 Mg/ha in mean 
carbon storage in all 15 sub-watersheds. 3.3. Timber volume in 2003 and 
urbanized 2060 Urbanization in the Lower Suwannee watershed produced a 
decrease of 2,300,000 m3 in the timber net volume in 2060, as a result of 
the loss of 22,600 ha of forest (Table 1); this despite an increase in 
timber volume of 1.7% used in the analysis (Prestemon S. Delphin et al. / 
Land Use Policy 54 (2016) 188-199 193 Fig. 2. Range of change in mean 
carbon storage (A) and water yield (B) between 2003-2060, by sub¬ 
watershed, in the Lower Suwannee River, Florida U.S. 194 S. Delphin et 
al. / Land Use Policy 54 (2016) 188-199 Fig. 3. Range of change in mean 
total carbon (A) and mean water yield (B) between 2003 and 2060, by sub¬ 
watershed, in the Pensacola Bay watershed, Florida U.S. S. Delphin et al. 
/ Land Use Policy 54 (2016) 188-199 195 Table 1 Total carbon storage 
(1000 Mg), mean carbon storage (Mg/ha), timber net volume (1000 m3), mean 
timber net volume (m3/ha), water yield (103 m3) and mean water yield (mm) 
for the Lower Suwannee River (LS) and Pensacola Bay (PB) watersheds, 
Florida U.S. Watershed-year Total carbon Mean carbon Timber volume Mean 
timber volume Water yield Mean water yield LS-2003 36,587 196 20,400 
109.5 93,193 1290 LS-2060 30,638 187a 18,100 110.6 96,695 1340a PB-2003 
7998 148 3700 68.6 104,631 1500 PB-2060 5867 137a 2900 68.3 108,441 1600a 
a Statistically significant differences at 95% confidence level, and Abt, 
2002a). However, the difference was not statistically significant at a 
95% confidence level. Timber volume in the sawlog portion, in 2060 was 
reduced by 1,940,000 m3 (Table 1). At the sub-watershed level, the 
largest decrease (17%) in mean timber net volume was for the Waccasassa 
Slough sub-watershed (number 49, Fig. 2); that had a decrease in forest 
area of 14% and increase in urban areas of 40%. In Pensacola Bay, there 
was a reduction of 11,500 ha of forest cover that resulted in a decrease 
of 800,900 m3 of timber net volume in year 2060 (Table 1). Sawlog timber 
volume was reduced by 22% in 2060 as a result of the forest area 
reduction. In 2060 there were no statistically significant decreases in 
the mean timber net volume at the sub-watershed level. However, 2 of 15 
analyzed sub-watersheds did show a non-statistically significant 
increase. These two sub-watersheds were not significantly affected by the 
projected urban areas. The largest difference (6 m3/ha) was in the 
Williams Creek-Oriole Beach Frontal sub-watershed (Id 11, Fig. 3); a 
decrease of 34% in forest area and an increase of 36% in urban areas. 

3.4. Water yield in 2003 and urbanized 2060 Mean annual water yield for 
2003 in the Lower Suwannee River watershed was 1290 mm and water yield 



volume was 12,900 m3/ha. However, in 2060, the mean simulated 
precipitation and PET were higher as were volume (13,400 m3/ha) and mean 
annual water yield (1340 mm);the latter being statistically greater. At 
the sub-watershed level, there were also increases in water yield, 
precipitation, and PET in 2060. The largest difference in mean annual 
water yield during the analysis period were in the Barnett Creek sub¬ 
watershed (number 1, Fig. 2) which also exhibited an increase of 152 
mm/year and the largest increase in precipitation in 2060; however no 
discernable changes in LULC were observed. In the 63 sub-watersheds there 
was statistically significant increased water yield per hectare in 2060. 
We note that uncertainties in future climate multi-decadal projections 
are currently large, and thus might mask an important effect of this land 
use change. There were statistically significant differences in the mean 
water yield between sub-watersheds in 2003 and 2060. In 2003,the annual 
water yieldinthe PensacolaBay was 1,520 mmandthe average water yield per 
hectare was 15,300 m3. However, in 2060 the assumed model precipitation 
and PET were higher as were annual water yield at 1,600 mm and water 
yield per hectare at 15,800 m3. At the sub-watershed level, the largest 
increase in simulated mean annual water yield was observed in the Bayou 
Cico-Bayou Garcon Frontal sub-watershed (number 15, Fig. 3) that 
registered the largest increase in precipitation in 2060 and there was a 
high positive correlation (r = 0.969) between the difference in annual 
water yield and difference in precipitation. There was also a decrease in 
the forest area and an increase in the urban areas. We assessed the 
InVEST model's water yield estimates in the Lower Suwannee River 
watershed to measure overall performance. Our assessment shows that there 
were differences between the average water volume at the Gopher gauging 
station and modeled average water yield. The average value from the model 
was 5,200,000 x 103 m3 while the average from the gauging station was 
6,700,000 x 103 m3. However, when values were normalized based on mean 
and standard deviation, trends between modeled and measured water yield 
curves were similar (Delphin et al., 2014). 3.5. Three-way interactions 
between carbon storage, timber volume, and water yield Although timber 
volume, carbon storage, and water yield provision levels are provided in 
Table 1 and Figs. 2 and 3, Fig. 4 shows the magnitude in interactions or 
trade-offs and synergies among these three services. In the 2060 
urbanized scenario, carbon storage and timber volume interactions were 
lower than in 2003; however, trade-offs with water yield were higher in 
2060 (Fig. 4). Overall there was a positive correlation between carbon 
storage and timber volume interactions (i.e., synergy), and a negative 
correlation with water yield interactions (i.e., trade-off). Because of 
forest loss as a direct result of urbanization, it can be surmised that 
both carbon storage and timber volume will decrease. However, water yield 
increased when urban areas increased, butthe correlation was weakly 
negative (Table 2). 4. Discussion 4.1. Land use and land cover map in 
2003 and 2060 We have used an integrated system of models and geospatial 
data to demonstrate a more robust approach to projecting the spatio- 
temporal interactions of ecosystem services in two different subtropical 
watersheds. We have also identified the areas in both watersheds that 
were most susceptible to ecosystem service drivers from projected urban 
development. Results show that integrating available geospatial models 
and forest inventory and monitoring data can also be used to project the 
ecosystem service tradeoffs that occur across space and time for 
different land use policy scenarios. Such knowledge on how land use 
change alters ecosystems services is fundamental for effective policy 
making and minimizing the negative effects on human well-being. But, we 
admit that predicting how future land use change might affect multiple 
ecosystem services is a difficult and complex task. And thus, no single 
existing predictive model is likely to be effective at addressing this 
question. Urbanization of existing forest cover types - not accounting 



for urban forest ecosystem service provision - directly affected carbon 
storage, timber volume, and water yield. Indeed urbanization is the 
dominant land use change that affects Southern U.S. forests (Alig et al., 
2004; Conner and Sheffield, 2006; Nowak et al., 2005; D N Wear and Greis, 
2002). Wear and Greis (2002) report that 4-9 million hectares of land are 
likely to be converted to urban land uses by 2060 in the region alone. 
Thus, it is likely that urbanization will convert many peri-urban forests 
by 2060, and result in significant decreases in the provision of carbon 
storage (Timilsina et al., 2011) and timber volume; but increases in 
water yield. Since LULC - as opposed to seasonal precipitation - input 
parameters are very sensitive in the InVEST model, one can surmise that 
urbanization is 196 S. Delphin et al. / Land Use Policy 54 (2016) 188-199 
Fig. 4. Ecosystem service interactions among carbon storage, timber 
volume, and water yield during 2003 - 2060 for the Lower Suwannee (LS; in 
black bars) and Pensacola Bay (PS; in grey bars) watersheds, Florida U.S. 
Table 2 Pearson Correlation coefficients between carbon storage, timber 
volume, water yield and percent (%) forest change during 2003-2060 at the 
sub-watershed level for the Lower Suwannee River (LS) and Pensacola Bay 
(PB) watersheds, Florida U.S. %Forest change Carbon Timber Water LS 
%Forest change 1 LS Carbon 0.98 1 LS Timber 0.99 0.98 1 LS Water -0.86 - 
0.93 -0.87 1 PB %Forest change 1 PB Carbon 0.99 1 PB Timber 0.97 0.97 1 
PB Water -0.14 -0.13 -0.20 1 more influential than precipitation in 
affecting water yield in this analysis (Sanchez-Canales et al., 2012) 

4.2. Carbon storage in 2003 and urbanized 2060 Mean carbon storage for 
the Lower Suwannee and Pensacola Bay watersheds were 196 and 148 Mg/ha, 
respectively, for the year 2003 and were within the range reported by 
Heath et al. (2011) for forest areas in the United States (162 Mg C/ha 

and a range of 150-192 Mg/ha). According to Timilsina et al. (2011) mean 

carbon stores for Florida forests were 165 Mg of carbon per ha in cold 
spots (i.e., areas of below average carbon storage values) and 177 Mg of 
carbon per ha in hot spots (i.e., areas of above average carbon storage 
values). As stated previously, carbon stored in urban areas was nottaken 
into account. Escobedo et al.(2009)presenteda comparisonof average tree 
carbon storage densities by urban forests in Florida that ranged 27.4 Mg 
per hectare in Pensacola and 30.8 Mg per hectare in Gainesville. Thus, 
using this study's 2003 average forest carbon storage values, the average 
carbon storage for Pensacola's urban forests (Escobedo et al., 2009), and 
adding this value to the mean carbon storage estimated in our study,the 
Pensacola Bay watershed will store about 175.4 Mg of carbon per ha 
representing about 16% of carbon stored by the watershed's urban forests. 
Using carbon storage values that include both urban and natural forests, 
total carbon storage in the urbanized 2060 scenario in Pensacola Bay 
watershed will decrease by 26% when compared to the 2003 total carbon 
storage. Thus, including urban forest carbon stores in our analysis; 
there would still be an overall decrease in total carbon storage in the 
Pensacola Bay watershed from 2003 to 2060. In Pensacola Bay, the forest 
area loss due to urbanization was 21% and findings corroborate other 
studies that show that land use change is an important driver affecting 
carbon storage and sequestration in ecosystems (Birdsey et al., 2006; 
Moore et al., 2002; Reyers et al., 2009), in addition to diseases, pests, 
hurricanes, and wildfires (McNulty, 2002) . 4.3. Timber volume in 2003 and 
urbanized 2060 The Lower Suwannee had 11% less timber volume in 2060 
compared to 2003; however, Pensacola Bay had a larger decrease (21%). As 
previously mentioned, the Pensacola Bay is more influenced by 
urbanization. Indeed Conner and Hartsell(2002) reportthat Florida has 
lost most of its forest areas due to urbanization. Hodges et al. S. 
Delphin et al. / Land Use Policy 54 (2016) 188-199 197 (1998) found that 
the forest area was reduced by urbanization by 25% compared to this 
study's 12% and 21% in the Lower Suwannee and Pensacola Bay watersheds 
respectively. 4.4. Water yield in 2003 and urbanized 2060 In both 



watersheds, water yield was higher in 2060 as were precipitation and PET. 
Based on CSIRO Mk3 and other studies, both precipitation and 
evapotranspiration should increase during this study's analysis period 
(Alvarez et al., 2001; Florida Oceans and Coastal Council, 2009; Furniss 
et al., 2010). Increase in water yield is primarily due to increase in 
precipitation (Dymond et al., 2012; Sun et al., 2005) and changes in 
LULC. Again, studies have documented lower water yield in forested areas 
(Bosch and Hewlett, 1982; Dymond et al., 2012; Hibbert, 1965).Inmany sub¬ 
watersheds in our study, the increase in water yield was the result of 
the increase in precipitation, decrease in potential evapotranspiration, 
and decrease in forest areas (Dymond et al., 2012) . The increase in water 
yield can however also have negative effects (i.e., ecosystem 
disservices), such as increased flood risk, nutrient transport, decreased 
water quality, and sedimentation of water bodies. Differences in observed 
water yield between the InVEST model estimates and measured gauging 
station data have been previously observed by Delphinet al. (2014) . These 
differences could be related to groundwater flow since according to Katz 
et al. (1997) watersheds in the study area are characterized by the 

presence of karst features that enable ground and surface-water 
interaction. The fact thatthe InVEST model does not accountfor these 
groundwater processes, should be taken into consideration where model 
estimates are important for management and policy actions, thus further 
calibration of the InVEST model or use of other models is warranted 
(Delphin et al., 2014) . 4.5. Three-way interactions between carbon 
storage, timber volume, and water yield Trade-offs among ecosystem 
services occur when the increase in the provision of one service can 
decrease the provision of another service (Bottalico et al., 2016; 

Cademus et al., 2014; Rodriguez et al., 2006). In general, increases in 
provisioning services will result in a decrease in regulating services 
(Dymond et al., 2012; RaudseppHearne et al., 2010) . Therefore, as also 
reported by Cademus et al. (2014), when forests are managed for increased 

carbon storage, water yield can be a trade-off as observed in our study. 
Also, if timber production is prioritized, the timber-carbon relationship 
will be negative; as timber harvesting will result in a decrease in 
biomass and carbon storage in live trees. Similarly, an increase in water 
yield due to forest loss can be considered either a win-lose (i.e., 
trade-off) or win-win (i.e., synergy) situation depending on the 
socioeconomic value placed by society to specific ecosystem services in a 
watershed (Kreye et al., 2014) . For example if there is a demand for 
irrigation water, more timber will result in less water yield (win-loss); 
However, decreased timber from logging results in more water yield that 
can result in either more water for irrigation or increased runoff, 
sedimentation, and flooding (Dymond et al., 2012) . 4.6. Limitations and 
future research There are limitations to our modeling approach. First, 
the 2003 FVLC map has been revised to create a more current statewide 
Cooperative Land Cover, version 2.3 (Florida Natural Areas Inventory, 

2012). But because of the time of release it could not be used for this 
analysis. Also, the watersheds will likely experience interactions 
andfeedbacks that arenot captured whenmodeling for each ecosystem 
service. For example, non-linear threshold effects or urban development 
patterns are difficult to predict and model as are changes in future 
economic and social conditions at both the local and regional levels (Yue 
et al., 2013) . Additionally, there are major uncertainties about regional 
precipitation patterns and temperature projections associated 
withanthropogenic climate change (Giorgi and Francisco, 2000) . These 
direct effects and additional potential climate feedbacks to the 
vegetation, such as changes in fire risk and primary productivity, add to 
the uncertainty in predicting future changes in ecosystem service 
provision. Our approach and findings however indicate several potential 
areas for future research. First, other direct and indirect drivers such 



as wildfires, drought, pest outbreaks, sea level rise, climate change or 
even increased energy and real estate prices could be modeled 
individually or coupled to determine the effects on ecosystem service 
bundles of interest. Second, determining what key ecosystem services - or 
disservices- and their supply level are in demand (or not) by 
beneficiaries according to relevant scales and contexts is important. 

This is particularly useful for determining when drivers will 
resultintrade-offs or synergies and what specific economic valuation 
methods to apply. Finally, approaches similar to ours could be integrated 
into applied geospatial tools and modeling scenarios to better visualize 
the effect of drivers on ecosystem service outcomes. These could also be 
used for spatially explicit analyses of the role of conservation or land 
use policies and activities (e.g., vegetation removal, urban development, 
land use ordinances) on long-term ecosystem service provision. 5. 
Conclusions Our findings highlighted the influence of a direct driver, 
urbanization, on trade-offs and synergies among ecosystem services during 
key periods of past, existing and projected increased urbanization rates. 
But, even though our findings show that there were - and will be - trade¬ 
offs among ecosystem services, some planning and management alternatives 
can be used for reducing these negative interactions. Thus, a key 
challenge when formulating policies for the provision of ecosystem 
services is identifying, defining and targeting optimal provision levels. 
As such, understanding the response of ecosystems to different drivers - 
using enhanced modeling alternatives and data- is key to finding these 
optimal levels and approaches such as ours can be used in identifying the 
areas that are the least "win-win" and most "lose-lose" and least 
susceptible to detrimental urbanization effects. An integrated system of 
models, including InVEST, regional land cover change and climate model 
alternatives, coupled with spatial data fusion from soil surveys, forest 
inventories, land cover classifications, and climate data sets can be 
used to project changes in ecosystem services associated with 
urbanization. This integrated approach is not exclusive to forested 
subtropical watersheds and can thus be used to better understand and 
identify the effects of drivers on relevant ecosystem structure (e.g., 
forest cover, soils, density, etc.) and subsequent functions (e.g., 
biodiversity, water regulation, primary productivity) and services when 
formulating policies in other geographic contexts as well. The results 
can be useful for land planners and managers in identifying specific 
conservation areas and potentially susceptible areas to forest land cover 
losses caused by urbanization. In doing so, the effects on ecosystem 
services can be used to assess differentland use policies, management 
practices, and planning alternatives that minimize these effects. We hope 
that studies such as ours can be used to not only better understand - and 
identify - the effects of drivers on relevant ecosystem services, but 
also to visualize the effects of policy decisions using geospatial tools 
and thus encourage the public and private sectors to participate in land 
management decisions. 198 S. Delphin et al. / Land Use Policy 54 (2016) 
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definitions (Woudenberg et al., 2010). Carbon pools FIA carbon data 
Aboveground biomassCarbon aboveground + carbon understory aboveground 
Belowground biomassCarbon belowground + carbon understory belowground 
Soil organic matter Carbon soil Dead biomass Carbon litter + carbon 
standing dead + carbon down dead 



